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I.  INTRODUCTION 

The  problem  of  the  excessive  computational  requirements  of  the 
Kalman  filter  has  been  studied  extensively  in  recent  years.  For  most 
practical  applications,  a  solution  of  this  problem  is  the  critical  fac¬ 
tor  in  deciding  whether  or  not  the  filter  can  be  implemented  in  an 
operational  system. 

The  problem  becomes  more  acute  in  the  case  of  a  delayed-state 
Kalman  filter  due  to  the  increased  complexity  of  the  filter  equations 
introduced  by  accounting  for  the  connection  of  the  measurement  with  the 
previous  state. 

One  approach  to  the  solution  of  the  computational  problem  has 
been  to  reduce  the  dimension  of  the  state  and  estimate  only  those  state 
variables  that  are  of  prime  interest;  see  Levy  [l4],  Schmidt  and 
Lukesh  [25],  Simon  and  Stuuberud  [27].  The  state  variables  often  ne¬ 
glected  in  this  reduced  filter  method  are  the  correlated  noise  processes 
in  the  system.  They  are,  in  effect,  modeled  as  uncorrelated  noise 
processes  when  this  method  is  used.  Although  optimal,  the  method  pro¬ 
posed  by  Simon  and  Stubberud  [27]  is  subject  to  the  restriction  that  the 
reduction  in  the  dimension  cf  the  state  is  limited  by  the  dimension  of 
the  measurement  vector. 

In  another  approach,  Aoki  and  Huddle  [3]  constructed  a  so-called 
minimal  order  observer,  utilizing  the  observer  theory  of  Luenberger  [l5j , 
which  was  composed  of  two  distinct  sections.  The  first  represented  a 
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dynamic  subsystem  of  the  observed  system  whose  output,  under  proper  con™ 
ditions,  was  a  linear  transformation  of  the  observed  system  state  vector. 
An  estimate  of  the  observed  system  state  vector  was  constructed  by  mini¬ 
mizing  the  elements  of  the  estimation  error  covariance  matrix.  The 
method  requires  proper  specification  of  several  matrices  for  the 
steady-state  condition  of  tue  filter. 

Brock  and  Schmidt  [4]  proposed  using  preset  or  precomputed  gains 
based  upon  the  high  and  low  frequency  response  of  the  system  in  an  iner¬ 
tial  navigation  application.  Sims  and  Melsa  [28]  investigated  the 
effect  of  using  constant  and  exponential  gain  functions,  which  could  be 
precomputed,  based  on  the  results  of  simulation  studies.  Precomputation 
of  gains  suffers  from  the  requirement  of  a  priori  knowledge  of  the  sys¬ 
tem's  behavior  in  order  to  predict  adequately  the  form  the  gain  should 
take. 

The  method  to  be  applied  here  to  the  delayed-state  filter  case 
is  due  to  Joseph  [12]  and  was  formalized  by  Pentecost  [22],  The  advan¬ 
tage  of  this  method  is  that  the  computational  requirements  are  signifi¬ 
cantly  reduced,  by  partitioning  of  the  system  into  lower  dimensional 
subsystems,  with  only  minimum  a  priori  knowledge  of  the  system's  behav¬ 
ior;  i.e.,  the  gains  are  computed  on-line. 

The  degradation  in  the  performance  of  a  continuous  Kalman  filter 
introduced  by  the  use  of  a  continuous  reduced-order  filter  has  been 
studied  by  Huddle  and  Wismer  [11]  Nishimura  [16]  performed  a  sensitiv¬ 
ity  analysis  on  the  continuous  Kalman  filter  performance  for  inaccuracies 
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in  the  specification  of  the  noise  covariance  matrices  and  of  the  initial 
state  covariance  matrix. 

Griffin  and  Sage  [9]  studied  the  continuous  and  discrete  Kalman 
filter  and  the  continuous  smoothing  equations  for  sensitivity  to 
modeling  and  noise  covariance  matrix  errors  analogous  to  our  method  for 
studying  the  performance  of  a  suboptimal  discrete-time  delayed-state 
filter  due  to  simplified  modeling.  Large  and  small  scale  sensitivity 
indices  were  defined  also,  which  are  not  useful  for  our  situation  as 
they  are  only  practical  for  a  small  number  of  parameter  changes. 

Heffes  [10]  studied  the  sensitivity  of  the  gain  computation  to 
modeling  errors  in  the  discrete  Kalman  filter  and  Sorensen  [25]  deter¬ 
mined  bounds  on  the  estimation  error  covariance  matrix  of  the  optimal 
Kalman  filter. 

Price  [24]  derived  the  recursive  equations  for  the  suboptimal 
estimation  error  covariance  matrix  of  the  discrete-time  Kalman  filter 
using  a  perturbed  model  by  a  method  equivalent  to  the  one  used  here 
for  the  delayed-state  filter.  He  then  determined  bounds  on  the  sub¬ 
optimal  estimation  error  covariance  matrix  and  conditions  on  the  mode*, 
changes  under  which  it  was  asymptotically  stable  in  the  large,  using 
Lyapunov  theory  techniques.  The  determination  of  such  bounds  and  con¬ 
ditions  for  stability  are  not  considered  here  for  the  delayed-state 
filter.  A  method  similar  to  Price's  could  be  used  for  such  a  purpose. 

We  here  derive  expressions  for  the  estimation  error  covariance 
matrix  for  suboptimizations  of  the  delayed-state  filter  based  on 


modeling  simplifications  and  alternate  gain  computations  which  place 
less  stringent  computational  requirements  on  the  system  computer  than 
does  the  optimal  filter. 

A  method  based  upon  the  work  of  Pentecost  [22]  is  presented  for 
the  reduction  of  computational  requirements  of  the  delayed-state  filter. 

A  performance  index  is  defined  for  use  as  a  relative  measure  of 
the  performance  of  the  various  suboptimizations  compared  to  the  optimal 
filter.  The  results  of  a  simulation  of  an  integrated  inertial  navigation 
system  are  presented  as  an  application  of  the  theoretical  results. 
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II.  THE  DELAYED-STATE  KALMAN  FILTER 

As  previously  noted,  measurements  which  are  linearly  connected  to 
not  only  i.he  present  state  but  the  previous  state  as  well  require  a 
modification  of  the  usual  Kalman  equations.  This  was  first  done  by 
Brown  and  Hartman  [7'J.  The  purpose  of  th<s  chapter  is  to  summarize 
the  resulting  equations  and  introduce  the  notation  that  will  be  used 
in  the  sequel.  The  derivation  of  the  equations  will  be  presented  in 
the  appendix. 

The  random  process  to  be  estimated  is  assumed  to  be  a  linear  system 
satisfying  the  vector  differential  equation 

x(t)  =  A(t)x(t)  +  B  (t)w(t)  (2.1) 

where  x(t)  is  the  n-dimensional  state  vector  at  time  t,  A(t)  is  an  nxn 
time-varying  matrix  at  time  t,  w(t)  is  a  p-dimensional  white  noise 
input  vector,  and  B(t)  is  the  nxp  time-varying  matrix  which  connects 
the  input  to  the  state  at  time  t. 

Solution  of  (2.1)  results  in,  for  the  discrete- time  case, 

x(k+l)  =  $(k+l,k)x(k)  +  u(k)  (2.2) 

where  $(k+l,k)  is  the  nxn  state  transition  matrix  from  stage  k  to  stage 
k+1,  and  |u(k)|,  the  plant  noise,  is  a  white  noise  sequence  of  n-dimen¬ 
sional  vectors. 
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The  measurement  model  is  assumed  to  be  of  the  form 


y  (k+1 )  =  M(k+l)x(k+l)  +  N(k+l)x(k)  +  v(k+l) 


(2.3) 


where  y(k+l)  is  an  m-dimensional  measurement  vector,  M(k+1)  is  an  mxn 
time-varying  matrix  which  linearly  connects  the  present  state  to  the 
measurement,  N(k+1)  is  an  mxn  time-varying  matrix  which  linearly  connects 
the  previous  state  to  the  measurement,  and  -|v(k+l)^  is  an  m-dimensional 
white  noise  sequence,  th*>  measurement  noise. 

The  following  statistical  properties  of  the  system  (2.2),  (2.3)  are 
assumed: 

1)  |x(k),k  =  0,1,...  j.  and  |y(j),j  =  l,2,...|are  gaussian  with 


zero 


means . 


2)  |u(k),k  *  0,1,  ...|  and  |v(j),j  =  l,2,...|are  independent  gaus¬ 
sian  white  noise  with  zero  mean. 

3)  E'|x(j)u*  (k)|  =  0  tor  all  k  2  j,  j  =  0,1,...,  where  E^*j-  indi¬ 
cates  the  expectation  operator  and  the  prime  denotes  the  trans¬ 
pose. 

4)  Ew(j)u'  (k)>  -  0  for  all  k  £  j  +  1,  j  =  0,1,... 

5)  Eix(j)v' (k)i  =  0  for  all  j,k,  j  =  0,1,...,  k  =  1.2,... 

6)  E(y(j)v' (k)|  =  0  for  all  k  >  j  +  1,  j  =  1,2,... 

7)  E<v(j)v'(k)|  =  V(k)6^k,  j,k  =  1,2,...,  where  6^  is  the  Kronecker 

delta  function  and  V(k)  is  a  positive  semidefinite  mxm  matrix. 

8;  e|u(j)u' (k)|  =  H(k)5^k,  j,k  =  0,1,...,  where  H(k)  is  a  positive 
semidefinite  nxn  matrix. 
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The  recursive  equations  for  estimating  the  state  of  the  model  (2.2)- 
(2.3)  at  stage  k  can  be  summarized  from  Brown  and  Hartman  as 

Q(k)  =  M(k)P(k|k-l)M’ (k)  +  V(k)  +  N (k)P(k-l (k-l)N'  (k) 

+  M(k)  $(.:,k,-l)P(k-l  |k-l)N'  (k) 

+  [M(k)$(k,k-l)P(k-l|k-l)N' (k)l’  (2.4) 

K(k)  =  [P(k|k-1)M' (k)  +  §(k,k-l)P(k-l|k-l)N' (k) 3q_1  (k)  (2.b) 

P(k|k)  =  P(k jk-1)  -  K(k)Q(k)K' (k)  (2.6) 

x(k|k)  =  x(k|k-l)  +K(k)[y(k)  -  M(k)x(kjk-1)  -  N (k)x(k-l  |k-l)  ] 

(2.7) 

x(k+l|k)  =  $(k+l,k)x(k|k)  (2.8) 

P(k+l|k)  =  i(k+l,k)P(klk)§' (k+J.,k)  +  H(k)  (2.9) 

where  x(k|k)  is  the  optimal  estimate  of  the  stage  at  stage  k  given 
measurements  through  stage  k  and  x(k+l|k)  is  the  optimal  estimate  of 
the  state  at  stage  k+1  given  measurements  through  only  stage  k.  K(k) 
is  the  optimal  gain  matrix  at  stage  k.  P(k|k)  is  the  covariance  matrix 
of  the  estj.mation  error  x(k|k)  =  [x(k)  -  x(k|k)]  and  P(k+l|k)  is  the 
covariance  matrix  of  the  estimation  error  x(k+l|k)  =  Lx(k+1)  -  x(k+l|k)]. 
Q(k)  is  the  covariance  matrix  of  tb*-  measurement  estimation  error 
y (k |k-l)  =  [y(k)  -  y(k|k-l)]. 

These  recursive  equations  uniquely  specify  a  solution  at  any  time 
k  given  the  initial  conditions,  P(0|0),  and  x(0|0).  A  block  diagram 
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indicating  the  structure  of  the  delayed-state  filter  is  shown  in  Fig.  1, 
where  the  double  lines  denote  vector  signals  and  the  blocks  denote  ma¬ 
trix  operations. 
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III.  PERFORMANCE  ANALYSIS  OF  THE  DELAYED-STATE  FILTER 


When  an  analysis  on  a  dynamical  system  is  performed  for  which  the 
estimate  of  the  state  variables  will  ultimately  be  required,  the  ques¬ 
tion  arises  as  to  how  complete  the  model  must  be  to  adequately  predict 
the  behavior  of  the  system  and  yet  result  in  computationally  efficient 
estimation  algorithms.  As  an  aid  in  performing  such  an  analysis,  it  is 
useful  to  have  an  algorithm  for  determining  the  degradation  in  perfor¬ 
mance  incurred  in  using  a  simplified  model  of  the  system  dynamics. 

Also,  since  the  gain  computation  in  the  Kalman  filter  places  the 
most  severe  requirements  on  the  computational  abilities  of  the  computer 
in  an  operational  system,  it  is  sometimes  desirable  to  use  suboptimal 
gains  which  are  easier  to  compute.  Thus,  it  is  desirable  to  have  an 
algorithm  to  evaluate  the  degradation  of  the  estimates  as  compared  to 
those  estimates  achieved  with  the  optimal  gain  computation. 

The  purpose  of  this  chapter  is  to  derive  algorithms  which  accom¬ 
plish  both  of  the  above  purposes.  The  method  of  attack  for  the  simpli¬ 
fied  model  case  is  analogous  to,  but  different  than,  the  method  used 
by  Price  [24j.  For  use  as  a  figure  of  merit  for  the  various  models  and 
suboptimal  filters  that  may  be  used,  a  performance  index  will  be  defined. 
The  effect  of  modeling  differences  in  the  noise  covariances  are  not  in¬ 
cluded  here  but  could  easily  be  included  by  suitable  modification  of  the 


equations  that  follow. 
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A.  Performance  Analysis  of  a  Suboptimal  Filter  Using 

A  Simplified  Model 

Suppose  that  a  change  is  made  in  the  modeling  of  the  system  dynam¬ 
ics  for  some  reason,  such  as  reducing  small  terms  to  zero  or  making 
simplifying  assumptions  which  result  in  changes  in  some  terms  or  reducing 
the  number  of  state  variables.  The  result  is  a  perturbation  in  the 
A (t)  matrix  of  the  system  model  (2.1)-(2.3)  which  introduces  a  perturba¬ 
tion  in  the  filter  equations,  (2.4)- (2. 9).  If  the  original  model  is 
assumed  to  be  the  most  complete  and  correct  one  possible,  it  will  re¬ 
sult  in  optimal  estimates  by  the  filter.  Thus,  if  the  model  is 
changed  or  simplified,  suboptimal  estimates  will  result.  We  will  call 
this  model  the  design  model. 

Let  the  original  system  model  be  described  by  (from  chapter  II) 

x(k+l)  =  $(k+l,k)x(k)  +  u(k) 

(3.1) 

y  (k+1)  =  M(k+l)x(k+l)  +  N(k+l)x(k)  +  v(k+l) 
and  the  design  system  model  by 

x,  (k+1)  =  $  (k+l,k)x  (k)  +  u(k) 

d  d  d  0.2) 

y,  (k+1)  =  M(k+l)x ,  (k+1)  +  N(k+l)x,  (k)  +  v(k+l) 
aad 

where  ^  =  $  +  6§,  The  perturbation  in  the  original  system  transition 
matrix  61  results  from  computing  the  transition  matrix  using  the 
design  system  model,  A(t)  +  6a (t).  In  general,  there  is  no  simple  rela¬ 
tionship  between  6a (t)  and  6§(k+l,k)  and  one  must  compute  both 
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transition  matrices  to  get  the  value  of  S$(k+l,k)  at  any  stage  k. 

By  using  the  design  system  model  in  the  optimal  filter  equations, 

A 

we  obtain  an  optimal  estimate  x,  of  x,,  if  y,  were  the  available  mea- 

a  a  a 

surement.  The  estimation  error  covariance  matrix  for  estimating  x^ 
given  y^  is 

Pd  (k|k)  =  E-j [xd (k)  -  xd(k|k^xd(k)  -  xd(k|k)]'|*  (3.3) 

We  wish  to  use  the  gain  obtained  from  the  design  system  optimal 
filter  and  the  design  system  transition  matrix  to  obtain  a  (suboptimal) 
estimate  of  x.  The  estimator  is  of  the  form 

Xs  (k |k)  =xs(k|k-l)  +  Kd(k)[y(k)  -  (k |k-l)  3.  (3.4) 

The  only  difference  between  this  estimator  and  the  one  for  the  optimal 
filter  for  the  design  system  is  that  y  is  used  instead  of  yd  since  yd  is 
not  physically  available.  x.(k|k-l)  is  computed  as  $. (lc,k-l)x  (k-l|k-l). 

b  OS 

The  Pd  used  in  computing  Kd  is  not  the  estimation  error  covariance 
matrix  for  the  suboptimal  filter,  but  is  the  estimation  error  covariance 
matrix  for  estimating  xd.  The  estimation  error  covariance  matrix  for 
the  optimal  filter  for  the  original  system,  where  $  and  y  are  used  in 
the  filter  equations  of  chapter  II,  is 

P(k|k)  =  E-j[x(k)  -  x(k|k)][x(k)  -  x(k|k)]'|.  (3.5) 


The  estimation  error  covariance  matrix  for  the  suboptimal  filter,  using 
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§d  in  the  optimal  gain  computation  and  in  the  estimator  with  y,  is 

Ps(k|k)  =  E-|[x(k)  -  xg  (k|k)][x(k)  -  xs(k|k)]'|  (3.6) 

To  evaluate  any  degradation  in  performance,  we  wish  to  compare  P  and  P^. 
Thus  an  expression  for  Pg  is  required.  The  relationship  between  the 
filters  mentioned  above  is  represented  schematically  in  Fig.  2. 

In  order  to  evaluate  P  ,  we  note  the  following: 

xs(k+l|k)  =  §(k+l,k)xs(k|lc)  +  t$(k+l,k)x,  (k|k) ,  (3.7) 

x(k+l|k)  =  $(k+l,k)x  (k|k)  -  6§(k-H,k)x  (k|k)  +  u(k),  (3.8) 

s  s 

x  (k+l|k+l)  =  $(k+l,k)x  (k|k)  -  &§{k+l,k)x  (kjk)  +  u(k) 
s  s  s 

-  Kd(k+l)ys(k+l|k).  (3.9) 

where 

y  (k+l|k)  =  y  (k+1)  -  y(k+l|k) 

5 

=  [M(k+l)§(k+l,k)  +  N0c+l)]xs(k|k) 

-  M(k+l)&$(k+l,k)xp(k|k)  +  v(k+l).  (3.10) 

With  respect  to  the  estimation  error  covariance  matrix  derivation,  the 
chief  difference  between  the  optimal  filter  and  the  suboptimal  filter 
considered  here  is  that  x  satisfies  the  original  system  homogeneous 
state  equation  and  the  xg  does  not.  As  we  shall  see,  considerable 
complexity  results. 


original  system  model  and  the  optimal  filter  for  the  design  system  model 
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We  first  compute  the  a  priori  suboptimal  estimation  error  covariance 
matrix 

Ps(k+l|k)  =  £|xs(k+ijk)x^(k+llk)|. 

After  a  great  deal  of  manipulation  of  (3. 7)- (3. 10) ,  we  obtain 
Ps(k+l|k)  =  $(k+l,k)Ps(k|k)  *'(k+l,k) 

+  H(k)  +  &$(k+i,k)S(k)b$'(k+l,k) 

-i  (k+l,k)LR(k)  -  S(k)]6$'(k+l,k) 

-  j$(k+l,k)[R(k)  -  S(k)]6$'  (k+l,k)j.'.  (3.11) 

If  second-order  effects  are  neglected, 

Pg(k+l|k)  =  $(k+l,k)Ps(k|k)  f'(k+l,k)  +  H(k) 

-  *(k+l,k)[R(k)  -  S(k)]6#'(k+l,k) 

-  j$(k+l,k)[R(k)  -  S(k)Jd$,(k+l,k)|’.  (3.12) 

We  see  that  Ps(k+ljk)  is  of  the  same  form  as  the  a  priori  estimation 
error  covariance  matrix  for  the  optimal  filter  plus  some  correction 
terms  due  to  the  change  in  modeling. 

To  compute  Ps(k+l|k+l),  rewrite  xs(k+l|k+l)  as 

xs(k+l|k+l)  =  xs(k+l|k)  -  (k+l)yg  (k+1  |k)  (3.13) 
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and  ys(k+l|k)  as 

y  (k+l|k)  =  M(k+l)x  (k+l'k)  +  N  (k+l)x  (k  |k)  +  v(k+l).  (3.14) 

J  s  s  s 

Ps(k+l|k+l)  =  e|xs  (k+1  |k+l)x'  (k+1  jk+1 

Note  that  x^k+llk+l)  and  y"s  (k+1  [ k )  are  of  the  same  form  as  in  the  opti¬ 
mal  filter.  Thus,  the  first  few  steps  of  the  derivation  are  the  same 
as  in  the  optimal  filter.  The  similarity  stops  when  one  tries  to  use 
the  equation  for  the  gain  to  simplify  the  equation  to  the  form  of 
equation  (2.6).  The  problem  is  that  the  gain  is  computed  using  Pd>  not 

P  ,  which  of  course  is  due  to  the  suboptimization  of  the  filter, 
s’ 

After  some  routine  manipulation, 

P  (k+l|k+l)  =  P  (k+ljk)  -  K  (k+l)[M(k+l)P  (k+l|k) 

S  S  Q  “ 

+  N (k+l)Pg  (k |k)  $d (k+l,k)  ] 

Kd[M(k+l)Ps(k+l|k)  +  N(k+l)Ps(k|k)$^(k+l,k)]j-' 
+  Kd(k+l)jM(k+l)Ps(k+l  |k)M'  (k+1) 

+  M(kH)  $  .  (k+l,k)P  (k  |k)N'  (k+1) 

+  [M(k+l)$d (k+l,k)Pg (k|k)N*  (k+1)]' 

+  N(k+l)Ps(;:jk)N'  (k+1)  +  V (k+1  )|k^ (k+1).  (3.15) 

In  order  to  completely  specify  P  ,  recursive  relations  for  R(k)  and  S  (k) 

must  be  found  ,  where 

R(k)  =  e[x (k)x'  (k|k)3 
s 
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R(k)  =  Hk,k-l)[R(k-l)  -  S(k-l)]|l  +  ^i(k)[§(k,k-l) 

-  6§(k,k-l)]+  N(k)|'K^(k)j  -  $(k,k-l)S(k-l)Ud(k,k-l) 
+  KJ(k)M(k)6§(k>k-l)]'  +  $(k,k-l)[R*  (k-1) 

a 

-  S(k-l)][M(k)$(k,k-l)  +  N(k)]'K^(k)  .  (3.16) 


Neglecting  second-order  effects. 


S  (k)  =  E[xs(k|k)x^(k|k)] 

=  $d (k,k-l)S (k-1) (k,k-l) 

+  $d(k,k-l)R*  (k-l)[M(k)§(k,k-l)  +N(k)]’K^(k) 

-  §d(k,k-l)S(k-l)[H(k)§d(k,k-l)  +  N(k)]’K^(k) 

+  Kd (k)-^R’  (k- 1 ) [M (k) $ (k ,k- 1 )  +  N(k)] 

-  S(k-l)[M(k)$d(k,k-l)  +  N(k)]*},$d(k,k-1) 

[M(k)$(k,k-1)  +N(k)]P  (k-l|k-l)[M(k)$(k,k-l) 

s 

+  H(k)]*  +  [M(k)$(k,k-1)  +  N(k)][R(k-l) 

-  S(k-1)]6#' (k,k-l)M' (k)  -  |[M(k)$(k,k-l) 

+  N(k)][R(k-l)  -  S  (k-l)]6$*  (k,k-l)M'  (k)j' 


+  Kd(k) 


+  V(k) 


Kd  (k) . 


(3.17) 


In  order  for  the  above  equations  to  be  uniquely  specified,  initial  values 
of  R(k)  and  S (k)  are  required.  Since  x(0|0)  is  usually  chosen  as  zero, 
R(0)  and  S (0)  are  null  matrices. 
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From  the  implementation  viewpoint,  note  that  computation  of  Pg  is 
not  required  in  the  actual  operation  of  the  filter.  However,  evaluation 
of  Pg  can  be  used  to  study  the  performance  of  the  suboptimal  filter  com¬ 
pared  to  the  performance  of  the  optimal  filter.  In  some  cases,  when 
the  actual  state  trajectory  is  known  in  the  simulation  study,  evaluation 
of  Pg  may  not  be  necessary,  only  a  comparison  of  the  optimal  estimates 
and  the  suboptimal  estimates  with  the  actual  state  trajectory  may  be 
sufficient  to  study  the  relative  performance  of  the  filters. 

B.  Performance  Analysis  of  a  Suboptimal  Filter  Using 
A  Simplifed  Gain  Computation 

If  the  original  model  is  used,  but  a  suboptimal  gain  computation 
is  used  for  some  reason  or  another,  the  P(k|k;  given  by  equation  (2.6) 
is  not  the  covariance  matrix  of  the  a  posteriori  estimation  error.  The 
a  priori  covariance  matrix  of  the  estimation  error  is  still  given  by 
equation  (2.9)  and  the  covariance  matrix  of  the  measurement  estimation 

error  is  still  given  by  equation  (2.4). 

Let  Pg(k+l|k+l)  be  the  covariance  matrix  of  the  a  posteriori  esti¬ 
mation  error  resulting  from  using  a  suboptimal  gain,  Kg  (k+1)  =  K(k+1) 

+  &K(k+l),  in  the  optimal  filter  equations. 

P  (k+1 | k+1)  is  given  by 
s 

P  (k+1 1 k+1)  =  ECx  (k+l|k+l)'x'  (k+1  |k+l) 3 

S  S  S 

=  Cl  -  K  (k+l)M(k+l)]P  (k+l|k)[l  -  K  (k+l)M(k+l) ]' 
s  s  s 

-  K  (k+l)P  (k(k)t'  (k+l,k.)[l  -  K  (fc+l)M(k+l)3* 
s  s  s 


19 


-  (k+l)Ps(k|k)§'  (k+l,k)[l  -  Ks(k+l)M(k+l)],|' 

+  (k+l)[N (k+l)Pg  (k  }k)N ’  (k+1)  +  V(ic+1)  ]lT  (k+1) .  (3.18) 

Expanding  (3.18)  and  neglecting  second  order  effects, 

Pg  (k+1  lk+1 )  =  [i  -  K(k+l)M(k+l)]Ps(k+l|k)[l  -  K(k+l)M(k+l)]' 

+  K(k+l)Ps(kjk)#'  (k+l,k)[l  -  K(k+l)M(k+l)  ]' 

-  |K(k+l)?s(k|k)$' (k+l,k)[l  -  K(k-.-l)M(k+l)]'| ' 

+  K  (k+1 )  [N  (k+1 ) Pg  (k  jk)N 1  (k+1)  +  V(k+1)]K,  (k+1) 

-  6k (k+1 ) [M (k+1 ) P  (k+l|k)  +  §(k+l,k)P  (kJk)N'(k+l) 

s  s 

+  N(k+l)P  (k  1  k) $ '  (k+l,k)  ]  +  6K(k+l)Q  (k+1  )K' (k+1) 
s  s 

+  K(k+1)[Q  (k+1)  -  M(k+l)P  (k+l|k)M'  (k+l)]6K'  (k+1). 
s  s 

(3.19) 

As  indicated  above  Q  and  P  (k+1  |k)  are  given  by 

s  s 

Q  (k+1)  =  M(k+l)P  (k+l|k)M' (k+1)  +  V(k+1)  +  N(k+l)P  (k|k)N' (k+1) 
s  s  s 

+  M(k+l)$(k+l,k)Ps(kjk)N' (k+1) 

+  [M(k+l)$(k+l,k)Ps(k|k)N' (k+1)]'  (3.20) 

and 

Ps(k+l|k)  =  $(k+l,k)Ps(k|k)§' (k+l,k)  +  H(k).  (3.21) 

The  above  equations  can  be  used  to  study  how  suboptimal  a  particu¬ 
lar  filter  is  and  how  sensitive  the  quality  of  the  estimates  are  to 
various  methods  of  computing  the  gain  matrix. 


C.  Definition  of  a  Performance  Index 


By  using  the  results  of  section  A,  we  can  obtain  the  perturbation 
in  the  estimation  error  covariance  matrix  for  a  particular  design  model. 
Likewise,  using  the  results  obtained  in  the  last  section,  the  perturba¬ 
tion  in  the  estimation  error  covariance  matrix  can  be  computed  from  a 
knowledge  of  the  perturbation  in  the  optimal  gain  due  to  a  suboptimal 
gain  computation.  Knowledge  of  these  quantities  may  be  useful,  but 
their  magnitude  does  not  give  a  complete  picture  of  how  the  overall 
performance  of  the  suboptimal  filter  compares  to  the  optimal  filter. 

The  information  needed  is  contained  in  the  value  of  the  cost  functional 
for  the  two  cases. 

The  cost  functional  considered  here  is  the  mean  squared  estimation 
error.  Let  J  (k)  be  the  value  of  the  cost  functional  for  the  optimal 
filter  and  Jg (k)  be  the  value  of  the  cost  functional  for  the  suboptimal 
filter.  From  the  optimality  of  J,  we  know  that  Js(k)  ^  J(k)  for  all  k. 
Then 

J  (k)  -  tr  E|x(k|k)x'  (k 
-  tr  P(kjk) 

(k|k)j 

=  tr  P(k  Jk)  +  tr  &P(kik). 


J  (k)  =  tr  E 
s 
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Define  a  performance  index  p,(k)  as 

J  (k)  -  J(k) 

^(k)  =  . Tqo~  (3*22) 

tr  6p(k|k) 
tr  P(k|k) 

We  see  that  p  for  fixed  k  has  the  properties  that  it  is  zero  and  a 
minimum  when  the  suboptimal  filter  is  actually  optimal  and  is  a  mono- 
tonically  increasing  function  o:  tr  &P;  i.e.,  the  more  suboptimal  a 
filter  with  respect  to  the  given  cost  functional,  the  larger  the  value  of  p, . 

This  performance  index  could  be  used  to  evaluate  the  relative 
performance  of  a  suboptimal  filter  compared  to  the  optimal  one  and  to 
evaluate  the  relative  sensitivity  of  the  optimal  filter  to  modeling 
variations • 
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IV.  SUBOPTIMIZATION  OF  THE  DELAYED-STATE  FILTER 


In  this  chapter,  we  derive  a  suboptimal  delayed-stace  filter  based 
upon  a  method  originated  by  Joseph  Cl2]  and  Pentecost  [22]  The  basic 
idea  is  to  transform  the  original  system  into  a  collection  of  smaller 
subsystems.  Optimal  estimation  is  then  performed  on  each  of  the  sub¬ 
systems.  An  estimate  for  the  total  state  vector  is  then  reconstructed 
from  the  estimates  of  the  subsystem  state  vectors.  This  estimate  is 
suboptimal  since  the  subsystem  estimators  are  unable  to  make  full  use  of 
the  information  available  regarding  the  total  state,  because  such  infor¬ 
mation  is  distributed  among  the  subsystems.  In  practice,  the  estimates 
of  the  subsystem  state  vectors  are  not  ne  ded  to  reconstruct  the  total 
system  state  vector  estimate;  only  the  subsystem  gain  matrices  are 
needed. 

For  convenience,  we  repeat  the  total  system  model  equations  as 
given  in  chapter  II: 


x(k+l)  =  *(k+l,k)x(k)  +  u(k) 

y(k+l)  =  M(k+1  )x  (k--l)  +  N(k+l)x(k)  +  v(k+l). 


(4.1) 


The  previous  assumptions  about  this  model  made  in  Chapter  II  are 
assumed  here  also. 


fcH 

Let  n^  be  the  dimension  of  the  j  subsystem,  where  j^nj  s  n 


and 


.th 


r  is  the  number  of  subsystems.  Let  m.  be  the  dimension  of  the  j 

r  ^ 

subsystem  measurement  vector,  .Z^m.  ^  m.  Define  the  j  subsystem  state 
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th  'n 

vector  and  the  j  subsystem  measurement  vector  by 


(k)  =  C^x(k) 


(4.2) 


Tl.  (k)  =  D.y(k) 


where  the  linear  transformations  C.  and  D.  are  n.xn  and  m.xm  respec- 

3  3  J  3 

tively  and  have  maximal  rank.  We  choose  and  to  be  time-invariant 
to  avoid  any  more  complexity  in  the  filter  than  necessary. 

At  this  point,  we  could,  by  sheer  brute-force,  proceed  to  design 
the  optimal  subsystem  estimator  as  in  the  design  of  the  optimal  estima¬ 
tor  for  the  total  system.  However,  we  will  take  a  simpler  route  by 
manipulating  the  subsystem  equations  into  the  form  of  (4.1)  and  then 
make  use  of  the  results  already  available  in  the  appendix. 

From  (4.2),  we  approximate  x  and  y  by 


x(k)  =  c|5.(k) 


(4.3) 


y(k)  =  (k) 

where  the  superscript  plus  denotes  the  pseudo-inverse;  i.e., 

c+  =  c'.  (c.cl)"1. 

3  3  3  3 

Premultiplying  the  first  equation  of  (4.1)  by  and  the  second 
equation  by  and  substituting  in  (4.2)  and  (4.3),  we  obtain 


S.(k+1)  =  C  Ak+l,k)cT^(k)  +  CjU(k) 


(4.4) 


1L(k+l)  =  D.M(kH)C^.(k+l)  +  DjN(k+l)cj§j(k)  +  DjV(k+l) 


Denoting  M.,  N.,  u.,  and  v.  by 
3  3  3  3  J  3 


^  (k+1  ,k)  =  Cj$(k+l,k)C 


Mj (k+1)  =  D  M(k+l)C 

N  .  (k+1)  =  D.N(k+l)C+ 
J  J  J 

Uj  (k)  =  C ^11  (k) 
v j  (k+1 )  =  D^v(k+1) 


(4.5) 


we  see  that  the  subsystem  equations  are  of  the  same  form  as  (4.1). 


(k+1)  =  $j(k+l,k)§j(k)  +u  (k) 
lu  (k+1)  =  Mj(k+1)§  (k+1)  +  N  (k+l)§  (k)  +  Vj (k+1). 


(4.6) 


Now  that  the  subsystem  equations  have  been  put  into  the  proper  form, 

it  is  a  routine  matter  to  apply  the  results  of  Brown  and  Hartman  [7]  to 

til 

each  subsystem.  Fcr  the  j  subsystem,  the  optim?l  delayed-state  filter 
equations  are 

Qj  (k)  =  M  (k)P  (k|k-l)MJ (k)  +  Vj  (k)  +  N  (k)Pj(k-l|k-l)N^(k) 


+  M  (k) I  (k,k-l)P^ (k-1 |k-l)N* (k) 

+  [fk (k)$ ^ (k,k-l)P^ (k-1 |k-l)Nl (k) ]' 


(4.7) 


K j  (k)  =  [Pj  (k  Ik-HMj  (k)  +  *  (k,k-l)P  (k-1  |k-l )NJ  (k)  Jq"1  (k) 


(4.8) 


PJ  (k  |k)  =  P^kjk-l)  -  K  (k)Q  (k)K’  (k) 


(4.9) 
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S.(k|k)  =  I.  (k |k-l)  +  K.(k)[Tl.(k)  -  M.  (k)t.(klk-l)  -  N .  (k) t.  (k-1  jk-1) ] 

J  J  J  j  J  J  J  j 


§^(k+l|k)  =  l(k+l,k)l(k|k) 

P^(k+l|k)  =  $  (k+l,k)Pj(k|k)$  (k+l,k)  +  KL  (k) 


(4.10) 


(4.11) 


(4.12) 


where 


V  .  (k)  =  E  jv  .  (k)v'  (k)r  =  D.V(kjD'. 

J  (  J  J  •  J  J 

H .  (k)  =  E-fu  .  (k)u'.  (k) >  =  C.H(k)C'.. 
J  l  J  J  J  J  J 


We  assume  that  the  estimate  of  the  total  system  state  vector  can 
be  reconstructed  from  the  estimates  of  the  subsystem  state  vectors  by 
the  relations 


a  XT*  ^ 

j:(k |k-l)  -  >  F.5.(k|k-1) 

J  3 


(4.13) 


:|k)  =  Y  F  (k|k) 
i  ]  J 

j=l 


(4.14) 


where  the  F.  and  C.  must  satisfy 
J  J 


y  F.C.  = 

Z,  J  J 
j=l 


(4.15) 


in  order  that  the  total  stace  vector  is  reconstructed  from  the  subsystem 
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state  vectors. 


From  (4.14)  and  (4.10),  we  have 

r  r 

x(k|k)  =  S  F  .  '  .  (k  |k— 1)  +  y  F.K.(k)D.[y(k)  -  M(k)x(k|k-1) 
]  ]  is  ]  ]  J 

j=l  j=l 

-  N (k)x(k-l  1  k- 1)3 


=  x(k|k-l)  +  J  FjKj  (k)D.[y(k) 

j=l 

-  M(k)x(k |k-l)  -  N (k)x (k-1 |k— 1)3 


(4.16) 


A  comparison  of  (4.16)  and  (2.15)  shows  that  the  partitioned-state 
estimation  method  implies  that 


K(k)  =  y  F  .K.  (k)D .. 

J  J  J 


(4.17) 


It  is  now  obvious  from  equations  (4.16)  and  (4.17)  that  estimates 
of  the  subsystem  state  vectors  are  not  needed  in  order  to  estimate  the 
total  system  state  vector.  With  the  estimator  now  in  the  form  of  (4.16), 
the  constraint  (4.15)  is  no  longer  needed. 

Note  that  5^(0  |0)  =  D^x(ojo)  and 

p^olo)  =  e{[=  (0)  -  I  (0 10)][ 5  (0)  -  ^(0)0)]'} 


=  CjP(0|0)C^. 


(4.18) 


Thus,  the  suboptimal  filter  requires  the  same  initial  conditions  as  the 
optimal  filter.  A  schematic  diagram  of  the  suboptima1  filter  is  given 
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in  Fig.  3. 

Summarizing  the  suboptimal  filtering  equations,  for  j=l,...,r, 

x(kjk-l)  =  §(k,k-l)x(k-l |k-l)  (4.19) 

$.  (k,k-l)  =  $(k,k-])C*  (4.20) 

(k  |  k-1 )  =  $j(k,k-l)Pj(k-l|k-l)^(k,k-l)  +  H.(k-1)  (4.21) 

Qj(k)  =  M  (k)P  (k|k-l)Mj(k)  +  Vj (k)  +  (k)Pj (k-1 |k-l)Nj (k) 

+  Mj (k)  (k,k-l)P j (k-1 jk-l)Nj (k) 

+  [Mj(k)$.(k,k-l)Pj(k-l|k-l)N’(k)]'  (4.22) 

K.(k)  =  [P.(k|k-l)M;(k)  +  $ (k,k-l)P  (k-l|k-l)NUk)]Q"1(k) 

J  J  J  J  J  J  J 

(4.23) 

P.(k|k)  =  F^ (k|k-l)  -  Kj(k)Qj(k)K*(k)  (4.24) 

x(k|k)  =  x(kjk-l)  +  K(k)ly(k)  -  M(k)x(kjk-1) 

-  N(k)x(k-1 |k-l) ]  (4.25) 


where  K(k)  is  given  by  equation  (4.17). 

Equation  (4.24)  is  true  only  if  the  gain  used  is  the  optimal  gain 
for  the  j1"*1  subsystem  and  is  computed  exactly.  If  is  not  computed 
exactly  (truncation  errors  may  occur),  it  is  not  guaranteed  that 
Pj(k|k)  will  remain  positive  definite.  To  avoid  this  problem,  it  is 
desirable  to  have  an  alternate  way  of  computing  P^(k|k)  which  will 
assure  positive  definiteness.  The  required  form  is 
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P.(k|k)  =  j§.(k,k-l)  -  K  (k)[M  (k)#(k,k-l)  +  N  (k)]|p  (k-j|k-l)* 
|$j(k,k-J)  -  K.(k)CM.(k)$(k,k-l)  +Nj(k)]}' 


+  Kj(k)Vj(k)K^(k)  +  H(k-l). 


(4.26) 


By  using  the  above  equation  and  equation  (4.21),  the  a  priori  subsystem 
estimation  error,  covariance  matrix  can  be  eliminated  from  equations 
(4.22)  and  (4.23)  for  and  respectively.  is  then  given  by 

Q  (k)  =  M  (k)$  (k,k-l)P  (k-l|k-l)$'(k,k-l)M'(k) 

J  J  J  J  J  J 

+  M.(k)H.(k-l)M'.(k)  +  V.  (k)  +  N.(k)P.(k-l|k-l)N'(k) 

J  J  J  J  J  J  J 

% 

+  M  (k)§.(k,k-l)P.(k-llk-l)N'(k) 

J  J  J  J 


+  [M  (k)4 . (k,k-_)P  (k-1 |k-l)N’ (k) ]' 


(4.27) 


and  K_  is  given  by 


K .  (k)  =  j$  .(k,k-l)P.(k-l|k-l)[M.(k)  §.(k,k-l)  +N.(k)]’ 
J  \  J  J  J  J  J 

+  Hjfc-UM'OOJQ^OO  . 


(4.28) 


Notice  that  the  computations  in  equations  (4.21)- (4. 24)  involve 
matrices  of  order  less  than  mxm  and  nxn.  The  partitioning  of  the  sys¬ 
tem  into  subsystems  is  not  unique.  The  optimum  partitioning  for  a  par¬ 
ticular  problem  at  this  time  must  be  found  by  trial  and  error  and  phys¬ 
ical  intuition.  If  the  subsystem  state  vectors  are  not  cross-coupled 
with  any  other  of  the  subsystem  state  vectors,  either  through  the  plant 
or  measurement  equations  or  through  the  noise  processes,  then  this 
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scheme  would  be  optimal,  as  indicated  by  Aoki  [l]  and  Pentecost  [22], 

At  the  present  time,  there  is  no  optimal  method  for  choosing  the 
matrices.  Usually,  they  are  taken  to  be  C*. 

As  an  aid  in  evaluating  the  reduction  in  the  amount  of  required 
computation  achieved  by  the  use  of  the  above  results,  it  is  useful  to 
compute  the  number  of  multiplication  and  addition  operations  required 
in  one  step  by  the  filters. 

Assuming  the  inverse  is  computed  using  the  gaussian  elimination 
method,  the  number  of  multiplication  operations  Mq  and  the  number  of 
addition  operations  Aq  required  in  one  step  by  the  optimal  filter  equa¬ 
tions  given  in  chapter  II  are  given  by 


Mq  =  2n^  +  n2(l+7m)  +  nm(5+4m)  +  nf* 

A  =  2n^  +  7n2ra  -  3nm  +  3nm2  +  m(m2+2)  -  n. 
o 


(4.29) 


(4.30) 


The  number  of  multiplication  operations  Mg  and  the  number  of  addition 

operations  Ag  required  in  one  step  by  the  suboptimal  filter  equations 

(4.17),  (4. 19)- (4. 25)  are  given  by 

r 

M  =  n2+  5nm  +  [2n^  +  6n.m2  +  Sn^m.  +  n.nm  +  nul  (4.31) 

s  Z,  J  J  J  j  j  j  J 

j=l 


A  =  m(3n+l)  +  /  [2n"?  +  8n2m.  +  n.m.(4m.~7) 

j=l 

-  n2  +  n,m(m,-hn-l)  +  m.  +  nu]. 

J  j  j  3  3 


(4.32) 
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As  an  example,  consider  a  sixteen-dimensional  system  with  a  scalar 
measurement.  Three  subsystems  are  chosen,  two  are  three-dimensional  and 
one  is  twelve-dimensional.  The  scalar  measurement  is  used  in  all  three 


subsystems.  Then 

M  =  10385 

M  =  5615 

o 

s 

A  =  9971 

A  =  5011 

As  a  consequence  of  using  the  partitioned  suboptimal  filter,  a  56%  re¬ 
duction  in  the  number  of  multiplication  operations  is  achieved  and  a 
50%  reduction  in  the  number  of  addition  operations  required  in  one 
step.  The  substantial  reduction  in  the  multiplication  operations  is 
particularly  significant  since  multiplication  operations  are  much  more 
costly  in  terms  of  computer  execution  time  than  addition  operations. 

The  question  of  how  suboptimal  this  method  is  may  be  investigated 

in  the  following  manner.  Let  K  be  the  gain  matrix  of  the  optimal 

delayed-state  filter.  Then  the  perturbation  of  this  optimal  gain 

incurred  by  using  the  above  partitioning  scheme  is 

r 

6K(k)  =  V  F.K.(k)D.  -  K(k) .  (4.33) 

3  J  3 
3=1 

ok  thus  calculated  can  then  be  ..sed  in  the  performance  analysis  equations 
of  chapter  III  to  calculate  p,  the  performance  index,  p  is  a  relative 
measure  of  how  suboptimal  a  particular  partitioning  is  for  a  specific 
system  and  can  be  used  to  evaluate  several  partitionings  to  determine 
which  one  is  the  closest  to  being  optimal  (least  suboptimal). 
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V.  AN  EXAMPLE:  NNSS  INTEGRATED  INERTIAL/DOPPLER- 
SATELLITE  NAVIGATION  SYSTEM 


In  order  to  illustrate  how  the  preceding  results  may  be  applied, 
we  consider  an  inertial  navigation  system  augmented  by  the  Navy  Naviga¬ 
tion  Satellite  System  (NNSS),  also  referred  to  as  the  Transit  System. 

In  the  next  three  sections,  we  will  present  the  model  of  the  system  and 
derive  expressions  for  the  noise  covariance  matrices.  After  some  gener¬ 
al  comments  on  the  method  used  to  simulate  the  system,  we  will  present 
the  results  of  two  performance  analyses  for  modeling  changes  and  the  re¬ 
sults  on  the  performance  of  a  partitioned-state  suboptimal  filter. 

For  our  example,  we  consider  a  terrestial  (low-altitude)  vehicle 
using  a  strapped-down  inertial  system,  with  the  vertical  channel  being 
implemented  by  other  than  inertial  means.  A  geocentric  latitude-longi¬ 
tude  coordinate  system  is  used. 

Such  an  inertial  navigation  system  is  capable  of  providing  global 
navigational  information,  but  suffers  from  long  term  "drift"  errors.  To 
compensate  for  this  "drift,"  another  source  of  information  may  be  used  to 
periodically  correct  the  inertial  system.  That  is  the  role  of  the  NNSS 
here.  The  role  of  the  Kalman  filter  is  to  integrate  the  two  systems  in 
an  optimal  fashion.  A  schematic  representation  of  the  integrated  system 
is  shown  in  Fig.  A. 

Periodically  the  vehicle  receives  a  message  signal  from  one  of  the 
satellites.  The  message  signal  gives  the  satellite's  position  and  other 
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miscellaneous  information.  Due  to  the  relative  velocity  between  the 
satellite  and  the  vehicle,  a  doppler  shift  in  the  signal  frequency  oc¬ 
curs.  Using  the  satellite’s  position  and  a  measurement  of  the  accumu¬ 
lated  doppler  count  over  a  time  interval  (approximately  twenty  seconds), 
information  about  the  inertial  system's  errors  can  be  inferred.  Since 
the  measurement  process  is  noisy  and  the  errors  in  the  inertial  system 
are  inherently  random,  a  Kalman  filter  can  be  used  to  obtain  a  better 
estimate  of  the  inertial  system  errors  than  could  be  obtained  by  using 
the  raw  measurement. 

As  we  shall  see  in  Section  B,  a  delayed-state  filter  is  required 
because  the  measurement  vector  (a  scalar  in  this  case)  depends  upon  the 
number  of  doppler  counts  accumulated  over  the  interval  between  sample 
times.  The  fact  that  the  measurement  vector  is  a  scalar  makes  this  sys¬ 
tem  particularly  attractive  from  a  computational  point  of  view  since  the 
Q  matrix  is  a  scalar  and  taking  its  inverse  is  a  trivial  operation.  On 
the  other  hand,  a  delayed-state  filter,  which  places  greater  demands  on 
the  computational  ability  of  the  on-board  computer  than  the  usual  Kalman 
filter,  is  required. 


A.  The  Plant  Model 

The  error  model  for  the  inertial  system  may  be  described  by  two 
basic  equations.  Pitman  _  2 3 J .  The  two  equations  are 


+  at  X  .'  =  € 


(5.1) 
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6R+2tt)X6R+u)X6R+u>Xu)X6R  =  6a-<|»Xa  -  u>26R  /c 

o  tan  (5.2) 

where  all  the  above  variables  are  three-dimensional  vectors,  the  cross 
X  indicates  the  vector  cross  product,  and 

'1 1  =  4>  -  68, 

<P  =  platform  coordinate  frame  error  vector, 

68  =  computer  coordinate  frame  error  vector, 
u>  =  platform  angular  rate  vector  with  respect  to  an  inertial 
frame  of  reference, 

e  =  gyro  "drift"  rate  (bias)  error  vector, 

6R  =  radial  position  error  vector, 

6a  =  accelerometer  "bias"  error  vector, 
a  =  sensed  "acceleration"  vector  (including  both  inertial  and 
mass  attraction  acceleration) , 

6Rtan  =  tangential  component  of  R  to  the  earth, 

2 

u)^  =  ^/R,  gm  is  the  mass  attraction  acceleration  and  R  is  the 
nominal  magnitude  of  the  earth’s  radius  vector  plus  the 
nominal  altitude  of  the  vehicle. 

Equation  (5.1)  describes  the  "twenty- four  hour  dynamics"  in  terms 
of  the  difference  between  the  angular  error  in  the  orientation  of  the 
platform  frame  and  the  computer  frame.  Equation  (5.2)  describes  the 
position  error  propagation,  the  "eighty-four  minute"  or  "Shuler"  dynamics. 
The  third  component  equation  in  (5.2)  will  not  be  used  since  that  channel 
is  to  be  implemented  by  other  means,  such  as  by  an  altimeter. 
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Since  the  system  is  a  strapped-down  system,  the  gyro  drift  errors 
and  accelerometer  biases  given  in  the  platform  frame  must  be  transformed 
through  a  continuously  updated  direction  cosine  matrix  into  the  computer 
frame  to  be  used  in  equations  (5.1)  and  (5.2). 

The  random  process  driving  functions  in  the  plant  equations  (5.1) 
and  (5.2)  are  the  gyro  drift  rate  and  the  accelerometer  bias.  To  fit 
the  form  required  by  the  Kalman  filter,  these  processes  must  be  gaussian 
white  noise  with  zero  mean.  By  physical  observation,  this  is  obviously 
not.  true.  The  standard  method  for  circumventing  this  difficulty  will  be 
employed.  We  consider  the  components  of  £  and  5a  to  be  correlated  noise 
processes  driven  by  gaussian  white  noise  with  zero  mean. 

More  specifically,  we  assume  that  the  components  of  £  and  6a  are 
first  order  Markov  processes  of  the  form 

x  +  px  =  \/2ctV  f(t)  (5.3) 

where  f(t)  is  unity  white  noise.  Augmenting  the  state  with  six  addition¬ 
al  states  due  to  £  and  5a,  the  plant  model  is  now  in  the  form  required  by 
the  Kalman  filter. 

In  addition,  three  more  state  variables  are  needed  to  complete  the 
model.  One  is  needed  for  the  vertical  position  error,  one  for  the  verti¬ 
cal  velocity  error,  and  one  for  the  doppler  count  bias  error. 

The  vertical  velocity  error  and  the  count  bias  error  will  be  modeled 
by  an  equation  of  the  form  of  (5.3).  The  model  is  one  which  can  be  made 
to  fit  a  variety  of  situations  and  implementations  by  adjustment  of  the 
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parameters  a  and  (3.  For  very  small  8,  the  process  approaches  a  true  bias 
and  for  very  large  p  the  process  approaches  white  noise.  The  amplitude 
is  adjusted  by  a. 

'Hie  vertical  position  error  is  taker  to  be  the  integral  of  the  ver¬ 
tical  velocity  error  and  is  thus  a  non3tationary  process. 

The  coordinate  frames  to  be  used  are  defined  as  follows: 


X,  Y,  Z 


x,  y,  z 


x’,  y’, 


=  earth-fixed  frame  with  X  through  the  north  pole 
and  Z  through  the  intersection  of  the  equator  and 
the  Greenwich  meridian, 

=  geocentric  navigation  frame  with  x  north,  y  west, 
and  z  up 

=  body-mounted  instrument-cluster  reference  frame 
with  x1  nose,  y1  left  wing,  and  z'  through  the 
roof. 


The  relationships  of  the  three  frames  are  shown  in  Fig.  5. 


For  convenience,  we  replace  6R  and  6R  in  equation  (’  /)  by  their 

x  y 

angular  equivalents  in  terms  of  60^  and  60^, 


and  define  state  variables  as  follows: 
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*9  ■  €x'/Bo 

*10-  €y-/u,o 


xn1  =  e  |  /(l) 

11  z*  o 


*12  °  ax’/R“o 


x,  _  =  6a  ,/Ru) 
13  y  •  o 


x, ,  =  6a  ,  /Ruj 
14  z'  o 


’  2 

x.  c  =  6R  /Ruj 
15  z  o 


x16  =  6N 


where  6N  is  the  doppler  count  bias  error.  Note  that  the  state  variables 
have  been  chosen  to  be  dimensionless  quantities. 

After  a  great  deal  of  algebraic  manipulation,  we  obtain  the  follow¬ 
ing  nonzero  elements  of  the  A  matrix  (Brown  [5]): 


al,2  =tBz 


a.  „  =  -  to 

1,3  y 


a.  =  U)  C  , 
1,9  o  xx’ 


a,  =  u)  C  , 
1,10  o  xy' 
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=  (jo  C  , 

o  xz' 


=  0)  C 


=  OU  C  . 
o  yy' 


=  u)  C  , 
o  yz' 


=  («  C  , 
o  zx' 


=  u)  C  . 
o  zy* 


=  u)  C 


a/  c  =  u> 
4,5  o 


z  1  1  2  2  2 

-  — —  - —  +  —  (ou  -hu  -uu  ) 
it  Rou  ou  x  y  o 
o  o  J 


ou  R  ,  ,  ,uj  R,  , 

y  1  A(  y  )  1 

d~  +  n - rr* — L  +  —  ou  U) 

Rou  Rou  At  ou  x  z 

O  O  0 
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u)  R  A(u)  R)  (i)  u) 

z  z  x  y 

Ra>  Ru)  At  +  u) 
o  o  o 


z  .  1  ,  2.  2  2. 

- - —  +  -  lU+DHl) 

Ru)  At  u)  v  y  z  o' 
o  o 


Au)  cu  u) 
y  x  z 


u>  At  a> 
o  o 


a,  =  (i)  C  , 

7,12  o  xx’ 


a,  ,  0  =  (jo  C  . 

7,13  o  xy' 


a_  ,.  =  u>  C  , 

7,14  o  xz' 


a_  -  2» 

7,15  y 


a8,15  "  ®o 


a9,9  ”  ’  32 


a10,l()  _  *  ^3 
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all,ll  =  '  P4 

a12, 12  '  " 

a13, 13  =  '  P6 
a14,14  =  '  P7 
a15, 15  =  ‘  P1 
a16,16  =  '  P8 


(8  X  8) 


(8  X  8) 
diagonal 
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The  state  variables  were  ordered  so  as  to  simplify  the  computation  of 
the  transition  matrix. 


B.  The  Measurement  Model 

In  this  section,  we  derive  the  measurement  model  for  the  integrated 
NNSS.  The  approach  is  the  same  as  in  Brown  [5]  and  Brown  and  Hagerman  [61. 

The  equation  for  the  ideal  count  observed  at  the  receiver  for  the 
time  interval  (t^  ,,  t^)  is  given  in  Stansell  [30], 

N(k)  =  AFAT  +  [ p(k)  -  p(k-l)]  ,  (5.4) 

where 


AF  =  fixed  offset  in  transmitter  and  local  oscillator 
frequencies  (32  kHz), 

AT  =  time  interval  between  timing  marks  referred  to  the 
satellite  time  base, 

>  =  wavelength  of  the  local  oscillator, 

Vj 

o  =  range  from  the  observer  to  the  satellite. 

A  number  of  sources  of  error  are  present  in  the  count  measurement 
process,  ^irst,  oscillators  are  not  so  stable  that  the  local  oscillator 
offset  can  oe  maintained  at  32  kHz  indefinitely,  so  some  unknown  bias 
effect  will  be  present.  Secondly,  refraction  due  to  the  atmosphere  af¬ 
fects  the  rf  signal  from  the  satellite.  This  is  particularly  true  when 
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the  satellite  passes  low  over  the  horizon.  Lastly,  jitter  in  the  count¬ 
ing  electronics  introduces  additional  random  error. 

These  sources  of  error  in  the  count  measurement  will  be  considered 
as  the  sum  of  a  correlated  process  and  an  uncorrelated  process.  The 
correlated  part  we  model  as  a  first  order  Harkov  process  and  designate 
as  state  variable  x^.  The  uncorrelated  part  we  assume  to  be  the  gaus- 
sian  white  noise  process  in  our  measurement  model.  Thus,  the  measured 
count  is  given  by 

N  (k)  =  AFAT  +  r~  [p(k)  -  p(k-l)]  -  SN(k)  -  v(k) .  (5.5) 

m  Ag 

If  the  inertial  system  were  without  error,  the  computed  count 
would  be  given  by  equation  (5.4),  But  since  the  inertial  system  is  in 
error,  the  range  computed  will  be  in  error  by  some  amount.  Thus,  the 
computed  count  is  given  by 

N  (k)  =  AFAT  +  [p(k)  +  6p(k)  -  p(k-l)  -  6p(k-l)]  (5.6) 

C  XG 

where  6p  is  the  error  in  the  range  due  to  the  position  error  in  the  in¬ 
ertial  system.  We  assume  that  the  position  of  the  satellite  is  known 
perfectly. 

Taking  the  difference  between  N  and  N  ,  we  obtain 

cm 

y(k)  =  N  (k)  -  N  (k)  =  t—  [  6p(k)  -  6p(k-l)1  +  -5N(k)  +  v(k)  . 

C  m 


(5.7) 
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To  obtain  the  form  required  by  the  delayed-state  filter,  we  must 
examine  the  terms.  The  expression  for  6p  is  given  in  Brown  and 
Hagerman  [7]  as 


5n  =  —  i  (R-R  C  )*R  +  (RR  C  )50 
O'  s  zz  z  s  yz  x 

s  s 


<RRS  cxZ  >«v  <5-8> 

s 


where  R  is  the  radial  distance  from  the  center  of  the  earth  to  the  vehi¬ 
cle,  R  is  the  radial  distance  from  the  center  of  the  earth  to  the  satel- 
’  s 

lite,  and  C  ,  C  ,  and  C  are  the  direction  cosines  between  the 
xz  yz  zz 

s  s  s 

navigation  axes  and  the  satellite  zg  axis  which  coincideswith  the  radial 
vector  between  the  center  of  the  earth  and  the  satellite. 

The  range  is  computed  from 

c  =  (R2+R2-2  RR  C  )*  .  (5.9) 

s  s  zz 

s 

Noting  that  69  =  x.  ,  66  =  x, ,  6R  /’R  =  x„,  and  6N  =  x,,,  the  mea- 

x  4  y  6’  z  8  16’ 

surement  model  now  fits  the  required  form, 

y(k)  =  M(k)x(k)  +  N(k)x(k-1)  +  v(k) .  (5.10) 

where 


M(k)  =  [0  0  0  b(k)  0  c(k)  0  Ra(k)  0  0  0  0  0  0  0  l], 

N(k)  =[000  -b(k-l)  0  c(k-l)  0  -Ra(k-l)  0  0  0  0  0  0  0  0], 


a(k) 


R 

s 


b(k') 


G 
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c(k)  =  -  [  RR  C  ]  . 

>„p  s  xz 
G  s 

The  model  for  the  system  is  now  complete.  We  turn  now  to  the  de¬ 
rivation  of  the  noise  covariance  matrices  to  be  used  in  the  filter. 


C.  The  Noise  Covariance  Matrices 


The  derivation  of  the  nonzero  terms  in  the  H  matrix  is  a  straight 
forward  exercise  in  classical  random  process  theory.  The  value  of  the 
V  matrix  (a  scalar)  is  based  on  physical  intuition  and  simulation  results. 
Brown  [5]. 

By  inspection,  the  first  three  diagonal  terms  are  white  noise  con¬ 
stants  due  to  gyro  drift. 


hl,l  =U1 


h2,2  U2 


h3,3  U3 


Similarly,  the  fifth  and  seventh  diagonal  terms  are  white  noise 
constants  due  to  accelerometer  bias. 


h5,5  =  U5 


h7,7  =  U7 
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Tne 


state  equation  for  the  first  gyro  drift  component,  x^,  is 


*9(t)  =  -  32*9(t)  +  V2ap2  f(t)  . 


(5.11) 


Solving  equation  (5.11)  over  the  time  interval  (t^,  t^^) 


x^(k+l)  =  expC-S^At)  Xg(k)  +  J  V2 o^32  exp(-32T)  f(t-T)  dT  (5.12) 


where  At  =  t^+^  -  t^.  The  h^  ^(k)  term  is  the  covariance  of  the  response 
of  x^(k+l)  due  to  f(t)  between  t^  and  i.e.. 


h9,9(k) 


At  At 

=  E  if  V2ct202  exp(-32T)  f(t-T)  dT  f  exp(-g,s)  f(t-s) 

Jf\  -J  n 


it  it 


1 1 


cxp[-  $2(t+s)]e[  f  (t-T)  f(t-s)l  dTd.s 


At  At 


*2S2 


I  I  exp[-  02(t+s)]6 (T-s)  dTds 
^  O  ^  Q 


o2  [1  -  exp(-282At)] 


ds] 


where  *  is  used  here  to  denote  the  Dirac  delta  function. 
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Similarly: 

hiO  10(k)  =  CT3  "  exp(-2  P3At)] 
hil  il(k)  =  [1  -  exp(-2  g^At)] 

hl2  12^  =  °5  L1  *  exp(-2  ^At)”? 

h13  13 ^  =  ae  t1  _  exP<*2  36^fc)3 

h14  14  ^  =  a7  -1  ‘  exP(-2  P?At)] 

h15  15 (k)  =  [l  -  exp (-2  ^At)] 

h16  16^  =  a2  [1  -  exp(-2  SgAt)]  . 
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The  only  cross-covariance  terms  which  are  present  are  hg  ^  and 


h15  y  •ire  equal. 


h8,15(k)  =  h15,8(k) 


At  At 


/I  2a^Q[l  -  exp(-8,T)]  expC-p^s)  6(t-s)  dTds 

Jq 


2a  at  / 

-3 — -  <[!  -  expC-g^t)]  -  %[l  -  exp(-2P1At)j 


Of  course,  values  of  the  a's  and  the  p's  must  be  determined.  In 
our  simulations,  we  shall  assume  values  which  we  presume  correctly  char¬ 
acterize  the  random  processes.  We  will  base  our  choice  of  values  on 
simulations  of  this  system. 


D.  General  Comments  on  the  Simulation  Study 

A  few  comments  are  in  order  concerning  the  values  of  the  statistical 
parameters  used  in  our  simulations  and  how  they  were  obtained. 

The  value  of  P(0j0)  is  based  upon  the  discussion  in  Brown  [5]  on 
establishing  initial  cjnditions  for  the  optimal  filter  for  the  system 
considered  here.  The  actual  numerical  values  are  the  result  of  a  trial 
and  error  simulation  study  of  the  performance  of  the  optimal  filter  using 
actual  flight  test  data.  The  optimal  filter  performed  "best"  when  the 
chosen  p  irameters  were  used.  Upon  that  basts,  the  selected  values  are 
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assumed  to  correctly  characterize  the  random  processes  involved  for  the 
purposes  of  this  study. 

The  values  in  the  initial  covariance  matrix  are  chosen  in  the  xyz 
navigation  coordinate  system  and  are  then  referred,  via  a  linear  trans¬ 
formation  into  the  body-mounted  coordinate  system,  in  which  the  actual 
initial  alignment  of  the  system  takes  place. 

In  the  xyz  coordinate  system,  the  initial  estimation  error  covar- 
iance  matrix  is  a  diagonal  matrix  P  (0 ] 0)  with 
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The  linear  transformation  relating  the  state  in  the  xyz  coordinate 
system  to  the  state  in  the  x'y'z*  coordinate  system  is  given  by 


r8  I  0  1 

T  =  |  ’  Tf  [  a"l  0 
0  fa  1  "r“l 


where  F  is  the  three  by  three  direction  cosine  matrix  relating  the  angu¬ 
lar  orientation  of. the  body-mounted  coordinate  system  with  respect  to  the 
navigation  coordinate  frame. 


C  |  C  ,  C  | 

x'x  x  y  x  z 


r  = 


C  •  C  |  C  | 

y'x  y’y  y  z 


.  c  ,  c  f  c  , 

L  z'x  z’y  z'z 


The  initial  estimation  error  covariance  matrix  to  be  used  in  the 
filters  is  given,  in  the  x’y’z’  coordinate  system,  by 

P(0|0)  =  TP*(0|0)T' . 


The  variances  and  inverse  time  constants  used  to  characterize  the 
random  processes  in  the  system  model  are  as  follows: 


2  , — .2 

accelerometer  biases:  a  =  400  sec 
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The  measurement  noise  covariance  matrix  V(k)  was  set  to  the  constant 
value  of  1600.  1 

To  compare  the  performance  of  the  different  filters,  true  error 
curves  have  been  plotted  with  the  estimates.  The  true  error  curves  were 
obtained  by  tracking  the  vehicle  from  ground  stations.  Due  to  the  pro¬ 
prietary  nature  of  this  data,  the  plots  have  been  normalized. 
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E.  Design  Model  Based  on  a  Simplified  Direction  Cosine  Matrix 

In  order  for  the  A  matrix  in  the  system  model  to  be  evaluated  at 
each  stage,  the  direction  cosine  matrix  relating  the  body-mounted  co¬ 
ordinate  system  must  be  updated  at  each  stage.  Under  the  assumption  of 
level  flight,  the  computation  of  this  matrix  becomes  considerably  sim¬ 
plified.  We  discuss  here  the  results  of  using  a  design  system  model 
based  on  this  assumption. 

Define  the  roll  angle  a  as  the  rotation  of  the  body-mounted  system 
about  the  x  axis,  the  pitch  angle  3  as  the  rotation  of  the  body-mounted 
system  about  the  -y  axis,  and  the  yaw  angle  y  as  the  rotation  of  the 
body-mounted  system  about  the  -z  axis.  The  direction  cosine  matrix  T  ^ 
relating  the  body-mounted  system  to  the  navigational  system  becomes 
(see,  e.g.,  Pitman  [23]) 


cacy 

easy  -  sas3cy 

-sasy  -  cas3cy  “ 

r-i 

-c3sy 

cacy  +  sas3sv 

-sacy 

s3 

sac  3 

cas3sy 

where  c  and  s  indicate  the  cosine  and  sine  of  the  indicated  angle,  re¬ 
spectively.; 

If  level  flight  is  assumed,  the  only  angular  displacement  which  oc¬ 
curs  between  the  body-mounted  system  and  the  navigational  system  is  in 


VuW 


The  resulting  simplified  direction  cosine  matrix  is 
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Using  ^  *  in  the  evaluation  of  the  A  matrix  results  in  a  design 
model  which  is  used  in  the  delayed-state  filter  as  in  chapter  3,  section 
A.  We  call  this  model  design  model  1. 

Figure  6  shows  the  normalized  estimates  of  the  latitude  error  of  the 
optimal  filter  and  the  suboptimal  filter  resulting  from  this  design  model 
as  a  function  of  the  flight  time  in  Greenwich  mean  time.  For  comparison, 
the  normalized  true  latitude  error  curve  is  shown  also.  Similarly,  the 
longitude  error  estimates  are  shown  in  Fig.  7. 

Upon  comparison  of  Fig.  6  and  7,  we  see  little  change  in  the  filter’s 
performance  is  introduced  by  the  assumption  of  level  flight.  This  result 
was  found  to  be  true  for  the  estimates  of  the  other  state  variables  also. 

The  design  model  resulting  from  this  simplification  is  by  itself 
not  much  simpler  than  the  original  model.  However,  this  simplification 
would  probably  be  made  in  conjunction  with  other  modifications  in  the 
original  model  such  as  those  to  be  discussed  in  the  next  section.  To 
assess  the  effect  of  each  change  in  the  model  on  the  filter  performance, 
the  changes  must  be  made  one  at  a  time  in  order  to  determine  which  one 
is  at  fault  if  poor  performance  results. 
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Design  model  1,  normalized  longitude  error 
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F.  Design  Model  Based  on  Simplified  Dynamics  in  th  .-.evel  Channels 

We  present  here  the  results  of  using  a  design  model  based  upon  sim¬ 
plification  of  the  dynamical  equations  for  the  error  in  he  two  level 
channels. 

The  basic  equation  used  to  derive  the  two  level  channel  equations 
for  the  original  model  is  equation  (5.2). 

6R  +  2u)  X  6R  +  u)  X  6R  +  w  X  (o>  x  6R)  =  6a  -  i|r  X  a  -  w  R_ 

'  o  tan 

A  number  of  simplifying  assumptions  on  this  equation  seem  reasonable 
due  to  the  nature  of  the  v  hide  and  typical  flights. 

If  u)  is  set  equal  to  zero  (Aw  <  <  ^  in  a  At  interval) ,  the  centripital 
acceleration  is  neglected,  and  R  is  considered  to  remain  constant  (AR  <  <  R 
in  a  At  interval),  then  equation  (5.2)  becomes 


6R  +  2u)  X  6R  =  6a  -  ^  X  a  -  a)  R. 

o  tan 


and  the  two  level  channel  equations  become 


xc  =  -  m  x,  -  u>  x,  +  2u>  , 
5  o  1  o  4  z 


x  -  (u  [C  ,x. „  +  C  ,x.,  +  C  ,x-.]  -  2u;  x._ 
7  o  yx’  12  yv  13  yz'  14  z  15 


X ,  =  -  U)  X  -  2a>  X-  -  u>  x,  +  u>  [c  (X,,  +  C  f x  +  C  ,xw]  -  2o>  x, 

7  o2  z5  06  o  xx'  12  xy'  13  xz'  14  y  15 


P 1 1 1  U  HP* 
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By  replacing  tl.«s  elirrc>:.:s  of  Che  fifth  and  seventh  rows  of  the  A 
matrix  of  the  original  uodsl  with  the  coefficients  cf  the  above.  r:b-‘.otion3 
'.ie  obtain  the  design  model  A  matrix.  Wo  call  the  resulting  model  i 

mo«i el  2. 

Ihn  results  of  using  design  model  2  in  the  filter  ere  ahewr,  i 
Fig.  8  and  Fig.  9.  There  is  little  significant  difference  in  the  esti¬ 
mates  of  both  the  latitude  and  longitude  error.  In  fact,  the  suboptimal 
estimates  appear  to  be  a  little  better  than  the  optimal  ones.  At  first 
thought,  this  would  be  cause  for  concern,  until  it  is  noted  that  these 
results  are  only  one  sample  out  of  a  possibly  large  ensemble  of  samples. 
The  only  conclusion  to  be  drawn  is  that  the  above  assumptions  do  not 
significantly  alter  the  performance  of  the  filter,  at  least  for  flight 
trajectories  similar  to  the  one  simulated. 

The  results  of  incorporating  the  simplifications  of  design  model  1 
and  design  model  2  in  one  model  are  shown  in  Figs.  10  and  11.  The  effect 
was  nearly  the  same  as  using  only  one  or  the  other  of  the  simplifications; 
in  fact,  the  performance  was  slightly  better  in  tha  longitude  channel. 

With  efficient  programming,  a  significant  reduction  can  be  made  in 
the  computational  effort  required  to  compute  the  transition  matrix  with¬ 
out  significant  degradation  in  the  filter's  performance  by  incorporating 
the  simplifications  of  design  models  1  and  2.  This  is  important  because 
the  transition  matrix  must  be  computed  at  the  end  of  each  20-second  in¬ 
terval  throughout  the  duration  of  a  flight. 
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Fig,  10,  Design  models  1  and  2  combined,  normalized  latitude  error 
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G.  Suboptimization  Using  the  Partitioning  Method 


Attempts  were  made  to  apply  the  partitioning  method  of  chapter  IV 
to  the  above  system  for  two  choices  of  partitions.  Although  achieving 
a  considerable  reduction  in  the  computational  requirements,  both  choices 
resul  ed  in  unsatisfactory  performance  of  the  filter.  Both  the  estimator 
and  the  estimation  error  covariance  matrix  equations  diverged  rapidly 
from  the  optimal  filter  results  by  several  orders  of  magnitude  after  the 
satellite  pass  period.  During  the  pass  period  when  measurement  data  is 
available,  transitory  oscillations  were  observed. 

In  the  first  partitioning  that  was  tried,  the  system  was  partitioned 
into  three  subsystems.  Subsystem  one  contained  x^,  x^,  and  x,.,  i.e., 
level  channel  one  and  the  coupled  psi  variable.  Subsystem  two  contained 
x2,  Xg,  and  i.e.,  level  channel  two  and  the  coupled  psi  variable. 

The  third  subsystem  contained  x^,  x^,  x^,  Xg,  .  .  .  ,  x^,  i.e,,  the 
twenty-four  hour  dynamics,  the  altitude  channel,  and  the  first  order 
Markov  instrument  noise  processes.  The  measurement  y  was  used  in  all 
three  subsystems. 

The  second  partitioning  tried  consisted  of  two  subsystems.  The  first 
subt^stem  contained  the  first  eight  state  variables  and  x^,.  and  x^g.  The 
second  subsystem  contained  x^,.  .  .  ,  x^,  the  instrument  noise  processes. 
The  measurement  y  was  used  in  both  subsystems.  This  choice  is  equivalent 
to  modeling  the  gyro  biases  and  accelerometer  biases  as  vrtiite  noise.  Due 
to  the  form  of  the  measurement  matrices  M  and  N,  the  second  subsystem 
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state  variables  are  estimated  as  the  initial  estimate  projected  through 
the  transition  matrix. 

The  divergence  problem  may  be  due  to  one  or  all  of  several  factors. 
First,  the  nature  of  the  form  of  measurement  for  this  system  may  not  be 
well  suited  to  suboptimization  by  this  method.  Information  about  sixteen 
state  variables  must  be  derived  from  a  scalar  measurement.  This  raises 


the  question  of  whether  or  not  the  observability  properties  of  the  sy.ttem 
place  a  constraint  on  the  choice  of  partitioning  or,  more  fundamentally, 
whether  or  not  the  partitioning  method  is  applicable  to  the  system  at  all. 
However,  both  the  observability  of  the  above  system  and  the  relationship 
of  observability  with  the  partitioning  method  remain  unexplored  problems. 

A  second  important  consideration  is  the  fact  that  only  two  parti¬ 
tionings  were  tried.  It  is  believed  that  the  successful  application  of 
the  method  to  this  system  is  highly  dependent  on  the  proper  choice  of 
partitioning. 

A  third  factor,  although  not  a  dominant  one,  to  consider  is  the 
error  due  to  the  digital  confutations.  In  any  study  of  this  sort, 
truncation  error  and  other  numerical  errors  will  occur  to  some  degree. 

The  simulation  of  this  system  is  particularly  vulnerable  since  many  of 
the  matrices  involved  are  ill-conditioned  and  a  great  number  of  matrix 
operations  are  required. 

As  a  final  consideration,  in  the  derivation  of  the  partitioning 


method,  the  pseudo-inverse  of  the  CL's  was  used  to  approximate  the  total 
system  state  in  order  to  de-couple  the  subsystem  model  from  the  total 
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system  model.  For  many  cases  this  may  be  a  very  poor  approximation. 

For  example,  the  first  two  subsystems  of  the  first  partitioning  are 
only  three-dimensional  and  the  total  system  is  sixteen-dimensional. 

Three  of  the  total  system  state  variables  are  approximated  with  the  sub¬ 
system  state  variables  and  the  remaining  are  approximated  as  zero  for 
all  time.  As  indicated  above,  it  may  be  possible  to  improve  this  approx¬ 
imation  by  a  different  choice  of  partitioning. 
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VI.  CONCLUSIONS 

Many  attempts  to  solve  the  excessive  computational  problem  of  the 
Kalman  filter  may  be  categorized  into  two  classes  of  suboptimizations. 
The  first  class  may  be  characterized  by  a  simplification  in  the  system 
model  that  is  used  in  the  filter  equations.  The  second  class  takes 
advantage  of  a  suboptimal  gain  computation  with  the  other  filter  equa¬ 
tions  remaining  the  same  as  for  the  optimal  filter.  In  chapter  III, 
we  presented  a  performance  analysis  for  both  of  the  above  classes  of 
suboptimisation  for  the  delayed-state  Kalman  filter.  The  derivation  of 
a  suboptimization  of  the  delayed-state  filter,  which  belongs  to  the 
second  class,  was  presented  in  chapter  IV. 

In  chapter  V,  the  results  of  a  simulation  study  of  the  performance 
of  an  integrated  inertial/doppler-satellite  navigation  system  were  pre¬ 
sented.  It  was  found  that  for  two  types  of  model  simplifications  the 
suboptimal  filter  performed  very  well  in  comparison  to  the  optimal  fil¬ 
ter.  On  the  other  hand,  attempts  to  apply  the  results  of  chapter  IV 
to  that  system  were  not  very  successful. 

Several  problems  remain  unsolved.  Further  experimentation  might 
result  in  a  partitioning  which  results  in  a  satisfactory  performance  by 
the  partitioned  suboptimal  filter  for  the  above  example. 

A  simulation  of  the  performance  analysis  equations  and  use  of  the 
performance  index  was  not  attempted.  However,  the  performance  index 
was  computed  for  the  partitioned  filter  and  confirmed  the  other  indica¬ 
tions  that  the  filter  was  not  performing  satisfactorily.  Concerted 
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effort  may  result  in  the  reduction  of  the  performance  analysis  equations 
for  modeling  changes  to  a  simpler  and  more  practical  computational  form. 

Optimization  of  the  partitioned  filter  by  the  choice  of  the 
matrices  with  respect  to  the  mean-squared  estimation  error  appears 
possible.  However,  this  would  result  in  r  additional  recursive  equations. 
Their  dimensionality  would  be  less  than  that  of  the  total  system  state 
and  the  amount  of  total  computation  required  should  still  be  less  than 
that  required  by  the  optimal  filter.  The  partitioned  filter  would  still 
be  suboptimal  with  respect  to  the  optimal  one,  but  should  have  improved 
performance  over  the  partitioned  filter  presented  in  chapter  TV. 
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VIII.  APPENDIX: 

DERIVATION  OF  A  KALMAN  FILTER  WITH  DELAYED-STATES  AS  OBSERVABLES 

We  present  here  a  derivation  of  t:he  recursive  equations  for  the 
optimal  delayed-state  Kalman  filter.  The  results  of  which  were  first 
obtained  through  the  work  of  Brown  and  Hartman  [7].  Our  method  is  not 
the  only  approach  that  could  be  used.  The  derivation  is  presented 
merely  to  make  more  rigorous  and  complete  the  results  in  the  preceding 
chapters. 

The  optimality  of  the  filter  is  a  consequence  of  a  fundamental 
theorem  in  estimation  theory  due  to  Sherman  [26]  which  we  will  state 
without  proof  (see  Deutsch  [8]).  We  will  then  state  a  more  specialized 
result  which  will  be  used  in  our  derivation. 

Let  x(k)  be  an  n-dimensional  random  variable  with  m-dimensional 
measurements  y(l) ,  .  .  .  ,  y(j).  The  conditional  probability  distri¬ 
bution  function  of  x(k)  given  y(l),  .  .  .  ,  y(k)  is 

F[5|y(i),  ....  y(j)3  =  p[x(k)£ §[ z(j)] 

where  P  is  a  suitable  prooability  measure  and  z(j)  is  an  mj -dimensional 

vector  whose  components  are  y,(l),  .  .  .  ,  y  (1),  .  .  .  ,  y, (j),  .  .  .  , 

i  m  i 

Vj)* 

We  wish  to  estimate  x(k)  based  only  on  the  measurements  y(l),  .  .  .  , 
y(j).  We  denote  this  estimate  by  x(k|j).  We  are  concerned  here  only  with 
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the  filter  problem  when  j  -  k.  We  have  the  smoothing  problem  when  j  >  k 
and  the  prediction  problem  when  j  <  k. 

We  denote  the  estimation  error  by 

*(k| j)  =  x(k)  -  x(k| j). 

We  would  like  x(k)j)  to  be  zero.  When  it  is  not,  we  assign  a  penalty  or 
loss  for  an  incorrect  estimate.  We  do  this  by  specifying  an  admissible 
loss  function  Lfx(k]j)]  which  has  the  following  properties: 


The  following  theorem,  which  we  also  state  without  proof  (see, 
e,g,,  Meditch  [18 1  or  Kalman  ri3  ] > ,  is  a  direct  consequence  of  Theorem  1. 


It  is  this  result  upon  which  the  Kalman  filter  equations  are  based, 
THEOREM  2. 


If  only  the  first  and  second  moments  of  the  stochastic  processes 
|x(k) ,k  =  0,  1,  .  .  .j  and  jy(i),  i  =  1,  2,  .  .  .  ,  4  are  known,  then 
the  optimal  estimate  for  all  admissible  loss  functions  is  the  linear 


estimate 


x(k|j)  =  Ejx(k)j  +  PxzP^  |z(j)  -  E[z(j)]| 

where  P^^  is  the  n  x  mj  cross-covariance  matrix  of  x(k)  and  z(j)  and 

P  is  the  covariance  matrix  of  z( j) . 
zz  VJ/ 

The  model  of  our  dynamical  system  was  given  in  chapter  II.  We 
repeat  it  here  for  convenience. 

x(k+l)  =  4  (k+l,k)x(k)  +  u(k)  (A.  1) 

y  (kf  1)  =  M(k+l)x(k+l)  +  N(k+l)x(k)  +  v(k+l)  (A.  2) 

The  assumptions  placed  upon  this  model in  chapter  II  are  required  here 
also. 

From  Theorem  1,  the  optimal  estimate  at  k+1  given  measurements 
y(l)  »  ....  y(k+l)  is 

x(k+l| k+1)  =  E-jx(k+l)  ]z(k+l)|  .  (A.3) 

For  any  gaussian  zero-mean  random  variables  x,  y,  z  (Papoulis  [21] 


or  Meditch  L 18  1) 
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where  z  =  z  -  E^zly},  and 

4iy>?}  =  E{*|y}  +  E^irj. 

Applying  these  formulas  to  (A. 3) ,  we  find  that 


(k+l|k+l)  =  E  jx (k+1 )  | z  (k )|  +  e|x  (k+1)  jy  (k+1  |k)} 
=  x  (k+1  jk)  +  E  jx  (k+1 )  jy  (k+1  j k )} 


(A.  4) 


where 


y(k+l|k)  =  y  (k+1 )  -  Ejy  (k+1)  |z  (k)|  . 

Letting  y(k+l|k)  be  the  optimal  estimate  of  the  measurement  y(k+l) 
given  measurements  through  time  k,  we  have  y(k+l|k)  =  y(k+l)  -  y(k+l|k). 
Expanding  y(k+ljk),  noting  that  v(k+l)  is  independent  of  |y(l),. . . ,y(k)}. 

y  (k+1  |k)  =  e}m (k+1  )x  (k+1 )  +  N(k+l)x(k)  +  v(k+l)  |z(k)J 

=  M(k+l)x(k+l|k)  +  N(k+l)x(k|k)  (A.  5) 

Noting  that  x(k+l|k)  =  $(k+l,k)x(k|k),  equation  (A. 4)  becomes 

x(k+l|k+l)  =  $(k+l,k)x(k|k)  +  E-jx(k+l)  |y  (k+1  |k)}  .  (A. 6) 

Since  x(k+l)  and  y(k+ljk)  are  zero  mean  and  gaussian,  we  have  as  a 
consequence  of  Theorem  2 

E  jx(k+l)  jy  (k+1  jk)}  =  P^P^Ck+ljk)  (A.  7) 

where  is  the  cross -covariance  matrix  of  x(k+l)  and  y(k+l|k)  and 
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„  is  the  covariance  matrix  of  the  measurement  estimation  error 

y? 

y(k+l|k);  i.e., 

=  E|x(k+l)y'  (k+l|k)| 

P__  =  E|y(k+l|k)y' (k+l|k)| 

Defining  K(k+1)  =  P^P^.,  we  have 

E|x(k+l)|y(k+l|k)|  -  K(k+1)  [y(k+l)  -  M(k+1) $(k+l,k)x(k|k) 

-  N(k+l)x(k|k)] 

=  K(k+l)[M(k+l)$(k+l,k)x(k)  +  N(k+l)x(k) 

+  vCk+1)  -  M(k+l)$(k+l,k)x(k|k) 

-  N(k+l)x(k|k)] 

=  K(k+l)[M(k+l)$(k+l,k)  +  N(k+l)]x(k|k) 

+  K(k+l)v(k+l).  (A.  8) 

Thus,  the  form  of  equation  (A. 6)  which  corresponds  to  equation  (2.7)  in 

chapter  II  is 

x(k+l|k+l)  =  x(k+l)  +  K(k+I)[y(k+1>  -  M(k+l)x(k+l|k) 

-  N(k-*-l)x(kjk)]  .  (A. 9) 

To  evaluate  K(k+1), 

P^~  =  E|x(k+l)y'  (k+l|k)| 

=  Ei[x(k+l|k)  -  x(k+l|k)][M(k+l)x(k+ljk) 

+  N(k+l)x(kjk)  +  v(k+l)]'^ 

=  P(k+1  |k)M'  (k+1)  +  $(k+l,k)P(klk)N’  (k+1)  (A.10) 

where  we  have  used  the  fact  that  v(k+l)  is  independent  of  x(k|k)  and 
P(k+l|k)  by  definition  is  E  jx(k+l|k)x(k+l jk)|. 
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Evaluating  P~~, 

yy 


=  E^y  (k+1  S  k)y'  (k+l|k)j 

=  E|[M(k+l)x(k+l|k)  +  N(k+?  )x(k|k)  +  v(k+l)  ][M(k+l)x(k+l  |k) 
+  N(k+l)x(k|k)  +  v(k+l)3'| 

=  M(k+l)P(k+l  |k)M!  (k+1)  +  V(k-*-l)  +  N(k+l)P(k |k)N' (k+1) 


+  M(k+l)§(k+l,k)P(k|k)N' (k+1) 

+  [M(k+l)§(k+l,Jc)P(k|k)N'  (k+1)]' 


(A.  U) 


Upon  comparing  (A. 11)  with  equation  (2.4)  we  see  that  P^  is  Q(k+1). 
Thus,  the  optimal  gain  matrix  K(k+1)  is  given  by 

K(k+1)  =  [P(k+l|k)M' (k+1)  +  $(k+l,k)P(kjk)N' (k+1) ]q_1 (k+1) 

(A.  12) 

where  Q(k+1)  is  given  by  (A.l’ ).  The  inverse  will  always  exist  since 
V(k+1)  is  assumed  positive  definite. 

P(k+l|k)  is  given  by 


P(k+1  |k)  =  E-|x(k+l  |k)x'  (k+l|k)j 


»  $(k+l,k)P(k|k)$”  (k+l,k)  +  H(k) 


(A. 13) 


Expanding  x(k+l|k+l),  we  obtain 


x(k+l  | k+1)  =  [I  -  K (k+1  )M (k+1 )  Jx (k+i  j k)  -  K(k+1)N  (k+l)x(k  |k) 
-  K(k+l)v(k+l) . 


Noting  that  E^x(k+1 Jk)x' (k |k)|  =  $ (k+l,k)P(k |k) , 
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P (k+1 j  k+1 )  =  E  jx  (k+1 1 k+1 )x' (k+1 1 k+1  >j 

=  P(k+l|k)  -  K(k+l)M(t +l)P(k+l|k) 

-  Ck  (k+l)M(k+l)P(k+l  |k)  ] ' 

-  K(k+l)N(k+l)P(k|k)  4*  (k+l,k) 

-  [K(k+l)N(k+l)P(k|k)$’ (k+l,k)]’ 

+  K(k+l)[M(k+l)P(k+l|k)M' (k+1) 

+  M(k+1)  §(k+l,k)P(k|k)N*  (k+1) 

+  |  M(k+1)$  (k+l,k)P(k jk)N’  (k+l)J-' 

+  N(k+l)P(k|k)N' (k+1)  +  V (k+1) ]K* (k+1)  (A.14) 

After  some  manipulation  of  (A. 14)  and  the  use  of  equation  (A. 13), 

P(k+1  |k+l)  =  P(k+l|k)  -  K(k+l)Q(k+l)Kl  (k+1).  (A. 13) 

Note  that  x(k|k)  is  a  gauss ian  zero  mean  process  and,  hence, 

P(0|0)  =  E-|x(0)x'  (0)|  since  x(ojo)  =  E^x(0)|^  =  0. 

This  completes  the  derivation  of  the  optimal  delayed-state  filter. 
The  data,  beginning  at  time  k  =  1,  are  processed  by  cycling  through 
the  recursive  equations  (A. 13),  (A. 11),  (A. 12),  (A. 13),  and  (A. 9)  with 
initial  conditions  P(0|o)  =  e|x(0)x' (0)\  and  x(0|0)  =  0. 
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